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Dissipation in noisy chemical networks: The role of deficiency 

M. Polettini,^'!^ A. Wachtel/ j^and M. Esposito^ 

Complex Systems and Statistieal Meehanies, Physies and Materials Seienee Researeh Unit, University of Luxembourg, 
162a avenue de la Faiencerie, L-1511 Luxembourg (G. D. Luxembourg) 

We study the effect of intrinsic noise on the thermodynamic balance of complex chemical networks subtending 
cellular metabolism and gene regulation. A topological network property called deficiency, known to determine 
the possibility of complex behavior such as multistability and oscillations, is shown to also characterize the 
entropic balance. In particular, when deficiency is zero the average stochastic dissipation rate equals that of 
the corresponding deterministic model, where correlations are disregarded. In fact, dissipation can be reduced 
by the effect of noise, as occurs in a toy model of metabolism that we employ to illustrate our findings. This 
phenomenon highlights that there is a close interplay between deficiency and the activation of new dissipative 
pathways at low molecule numbers. 


I. INTRODUCTION 

Today, advanced methods in genomics and 
metabolomics allow to reconstruct the chemical 
networks (CN) describing the metabolism of complex 
organisms^. These reconstructions are graphical 
repositories of thousands of pathways, metabolites, and 
their stoichiometry. Much like heat engines, metabolism 
operates thermodynamic cycles far from equilibrium 
that transform low chemical potential environmental 
resources into valuable products, at the expense of high 
chemical potential waste. Unlike the working substance 
of heat engines (e.g. steam), some metabolites, enzymes 
and cofactors might reach very low concentrations. 
At this level intrinsic noise, due to discreteness and 
randomness of molecular collisions, enters into pla}!^. 
Suppression of noise and control of correlations in 
the abundance of regulatory molecules is crucial for 
the correct functioning of metabolic networkJ^K^. A 
stochastic description of dynamics and thermodynamics 
based on jump processes in molecules’ populations is 
then required. 

In this direction, the growing field of Stochas¬ 
tic Thermodynamics created the basis for a com¬ 
plete and consistent characterization of irreversibility 
in small nonequilibrium systems subject to fluctua¬ 
tions. Dissipation is quantified by the rate at which 
entropy is produced (EPR) and eventually delivered 
to the environmenfP. The theory has been applied 
to general such as those involved in gene 

regulatiorP^, cellular computatiorf^, copolymerizatioiJ^, 
kinetic proofreading!^, chemical switched!^, and signal 
transductiorP^. On the other hand, there is a growing 
body of mathematical literature linking a CN’s topol¬ 
ogy to its dynamics, and still bearing no thermody¬ 
namic interpretation. In particular, it has been under¬ 
stood that a topological number called deficiency sub¬ 
tends the onset of complex behavior, such as bistability 
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and oscillation^!^^^, which are the mechanisms of chemi¬ 
cal switches and clock^^. When intrinsic noise is impor¬ 
tant, a crucial result by Anderson, Craciun and Kurtz 
(ACK]P^ relates the deficiency of the CN to steady sta¬ 
tistical properties of the chemical mixture. 

In this paper we merge stochastic thermodynamics and 
deficiency theory, via the ACK theorem. We compare the 
behavior of an arbitrary CN subject to intrinsic noise and 
that of the corresponding deterministic model without 
noise, which follows deterministic rate equations where 
correlations between species are neglected. In the limit 
of large particle numbers the deterministic dynamics de¬ 
scribes the mode, i.e. the most typical behavior of the 
system. The difference between the stochastic and the 
deterministic EPR in the two cases, here named corre¬ 
lation EPR (previously known as fluctuating EPR, to¬ 
day ambiguous), is known to vanish at steady states for 
linear CNs where only input/output and conformational 
changes of a molecule are allowed, and reaction velocities 
are linear-affine in the molecules’ population^. 

The main result in this paper is to extend this obser¬ 
vation to nonlinear CNs with null deficiency at steady 
states, and to linear networks at all times. We rely on 
the following formula for the steady correlation EPR as 
the weighted difference between the mean and the mode 
of the reaction velocity v, 

correlation EPR = (meanu — most probable u) G, (1) 


where G is the free-enthalpy increase. Hence the correla¬ 
tion EPR might be interpreted as a measure of a system’s 
“propensity to complexity”. 

The plan of the paper is as follows. In Sec. |IIB| we pro¬ 
vide a simple definition of deficiency with the aid of a toy 
modei of metaboiism. More generaiiy, under the assump¬ 
tion that the iaw of mass-action hoids and that the mix¬ 
ture is weii-stirred, we iiiustrate the dy nami cs and ther¬ 
modynamics o f CN s, in the stochastic ( II C[ ) and in the 
deterministic ( IID[ ) settings. We then derive the above 
formula, and by virtue of the ACK theorem (whose proof 
we briefly sketch in Appendix]^ we draw our main con¬ 
clusion that the correlation EPR vanishes for networks 
with zero deficiency. Our toy model will finally serve as 
a testing ground. We employ it to illustrate through Eigs. 
ii the predictions of the ACK theorem. Incidentally, 
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the model displays a non-positive correlation EPR, some¬ 
what contrary to the intuition that “large variability is 
likely to [... ] increase metabolic burden’®. We give an 
explanation of this phenomenon in terms of the topology 
of the state space where stochastic population dynam¬ 
ics occurs, showing that when deficiency is nonzero, for 
low molecule numbers certain irreversible closed reaction 
pathways are switched off. 


II. SETUP 
A. Notation 

As customary in CN studies, we employ a rather com¬ 
pressed notation. Letting X be the vector of chemical 
species, a CN is depicted by a set of stoichiometric equa¬ 
tions 

i^+p • X V i>-p ■ X (2) 

k—p 

where vectors iz+p and U-p contain, respectively, the 
numbers of molecules of each species being consumed and 
produced by reaction p, and a • 6 is the scalar product. 
The stoichiometric vector is defined as Vp := iz_p — iz+p, 
and it describes the net increase of species’ populations. 
The stoichiometric matrix is the matrix that has the stoi¬ 
chiometric vectors as columns, V = (Vp)p>o. We assume 
that all reactions are strictly reversible, that is, ^±p > 0 - 
In sums index p spans over reactions in both di¬ 
rections, unless otherwise specified. Analytic operations 
between vectors are performed component-wise and im¬ 
ply the scalar product, e.g. := 

a • Inh := In 6 ^. Boltzmann’s constant is set to 

unity. 


B. From metabolism to deficiency 

Roughly speaking, the deficiency of a CN is the num¬ 
ber of “hidden” closed pathways, or thermdynamic cy¬ 
cles. Let us make this more precise with a simple model 
inspired W metabolism. Emphasis is on the cycle struc¬ 
ture (see^ for a formal introduction). The model reads 

0 -4 N 

N + toE A (m + n)E + W (3) 

nE + W A 0, 

where 0 signifies the “environment” as a whole. The first 
reaction introduces nutrients N. The second processes 
the nutrients with the aid of m tokens of energy E to pro¬ 
duce more tokens of energy and waste W, and the third 
delivers waste and excess energy to the environment. 

When all three reactions in the above network are per¬ 
formed in a pathway, a thermodynamic cycle is com¬ 
pleted, restoring all concentrations in the system to their 


initial value at the expense of irreversibly dissipated free 
enthalpy (entropy production). Correspondingly, the sto¬ 
ichiometric matrix 


+1 -1 0 
V = I 0 + 1-1 

0 +n —n 


(4) 


admits c = (1,1,1)^ as a right-null vector, Vc = d^. 

The crucial step to understand deficiency is to intro¬ 
duce a symbolic representation of the network in terms 
of complexes^ which are aggregates of species appearing 
as either reactants or products in a reaction. In our case, 
the complexes are Yi = 0 ,Y 2 = N,Y 3 = N + mE, Y 4 = 
{m + n)E + W, Y 5 = W + nE. We then obtain a repre¬ 
sentation of the CN as a graph by drawing each reaction 
as an edge connecting vertices given by the complexes. 

Eor m = 0, we notice that Y 2 = Y 3 and Y 4 = Y 5 and 
that a representation of the above network in terms of 
complexes is a graph consisting of one cycle: 



(5) 


Its topology is fully described by its incidence matrix 


d = 


-1 0 + 1 \ 
+ 1-10 
0 +1 -1 / 


(6) 


which admits one right null vector. 

Eor m > 0 we obtain the representation 


Y5^-Yi^^Y2 


Y3^-Y4 , (7) 


with incidence matrix 


d = 


+1 

0 

0 

0 


0 

0 

-1 

+1 

0 


0 

0 

0 

- 1 / 


(8) 


This graph has no cycles; in fact its incidence matrix 
admits no right-null vectors. 

The deficiency (5 of a CN is the number of indepen¬ 
dent closed reaction pathways that cannot be visualized 
as independent cycles in the graphical representation in 
terms of complexes, and thus in some sense are “hidden”. 
In our example when m = 0 then (5 = 0, otherwise the 
system is deficient, ^ = 1. Notice that null deficiency 
occurs when the autocatalytic mechanism of reaction 2 
is not present. 

The general recipe to calculate the deficiency is: (i) 
write down the stoichiometric matrix V of the network; 
(ii) write down the incidence matrix d of the graph where 
the reactions are arrows and complexes of reactants dis¬ 
tinct vertices of the graph; (iii) then the deficiency is 

S = dim ker V — dim ker 5 > 0 (9) 
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where dim ker calculates the dimension of the null space. 
The deficiency is non-negative. In fact one can write 


V = 


dY , 


( 10 ) 


where the entry {dY/dX)ij quantifies the amount of 
species Xi in complex Yj. Since by Eq. (10) a right-null 
vector of d is necessarily a right-null vector of V, then 
S>0. 


It can easily be proven that the EPR is non-negative, 
embodying the second law of thermodynamics. The log¬ 
arithmic term measures the thermodynamic cost of re¬ 
action p for a given X, and it quantifies the degree by 
which detailed balance is broken. 


D. Deterministic EPR 


C. Average stochastic EPR 


The setup of Markovian population dynamics of chem¬ 
ical species is as follows. The number of molecules in the 
reactor performs a jump process on the discrete lattice 
orthant Zxq of populations that, starting from the ini¬ 
tial state Xo, are reachable by a finite number of reac¬ 
tion^ According to the law of mass-action, transition 
X ^ X + Vp is performed at rate 

- Up)\' 

The probability (or ensemble) Pt{X) that X molecules 
are present in the reactor at time t obeys the Chemical 
Master Equation pt = Lpt with generator 


Lft(X) = -^[u+p(X)ft(X) 

P 

-v.p{X + Vp)pt{X + Vp) 


( 12 ) 


Multiplying by, and summing over X, one obtains for the 
mean populations 

|(X)t = ^Vp(up(X))i (13) 

P 

where the average (• )t is taken with respect to Pt{X). 
The equation is not closed, as it involves higher moments 
on the right-hand side. 

Eor finite Zxq , it can be proven that any ensemble sup¬ 
ported on Zxo evolves towards a unique steady ensemble 
Poo such that Tpoo = 0- We assume that for unbounded 
Zxq conditions are met by which at all times pt{X oo) 
decays fast enough (e.g. exponentially) so that no proba¬ 
bility leak to infinity occurs, and that a steady ensemble 
exists. 

In this framework, the average EPR characterizing the 
CN’s dissipation is defined a^ 


The corresponding deterministic model is obtained by 
neglecting correlations and higher cumulants, i.e. by 
replacing {X^p)t where is a large vol¬ 

ume parameter that makes x a continuous variable with 
the interpretaton of a concentration; in the following 
we will set Q = 1 for notational clarity and only re¬ 
sume proper scalings when studying the model systems 
in Sec. imi Also, in the large volume limit the approxi¬ 
mation Vp{x) ~ kpX^p is made. Then Eq. (13) yields the 
rate equatioiP 


= ( 15 ) 

p 


Again, we are interested in steady behavior, when the 
right-hand side vanishes. Importantly, while the Chemi¬ 
cal Master Equation admits one unique steady ensemble, 
the corresponding deterministic dynamics might admit 
none or several locally stable fixed points Xoo and more 
complicated ^enomenology such as limit cycles and frac¬ 
tal attractor^. Deterministic multistability corresponds 
to the steady ensemble being multimodal. Notice that 
X cannot be interpreted as a mean, as for bistable sys¬ 
tems the mean might be far from both stable fixed points. 
Rather, in a scaling limit with the system size, random 
jump processes can be shown to typically behave deter¬ 
ministically, as rigorously detailed in Ref.^^. 

In this setting, the deterministic EPR is defined a^^ 




'^Vp{xt)\n 

P 


v+p{xt) 

'^—p{Xt) 


> 0 . 


(16) 


The connection to free-energy differences and other ther¬ 
modynamic potentials in a nonequilibrium setting is de¬ 
tailed in Ref.^. 


III. RESULTS 
A. Theoretical 


CTt := 


^(t)p(X)ln 

n ' 


v+p{x)Mx) \ 
v.p{X + Wp)pt{X + Wp)/, 


> 0 

(14) 


Eirst, we re-work the above expressions for the deter¬ 
ministic and stochastic EPRs to make them closer one to 
another. Introducing the thermodynamie forees 


Gp:=ln^, (17) 

K>—p 

^ That is, Zxq '•= + Vn, n G > 0}, some¬ 

times called the stoichiometric compatibility class, compatible that measure the kinetic imbalance of reactions, with a 
with Xq. manipulations we can bring the deterministic EPR 
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FIG. 1. We consider a class of toy models for metabolism 
0 ^ N, N + mE ^ (m + 2)E, 2E ^ 0, for varying m. In this 
figure we compare stochastic and deterministic time evolution 
of nutrient and energy molecules in the model corresponding 
to m = 0, that has deficiency (5 = 0 (on the left), and in 
model corresponding to m = 3, with deficiency (5=1 (on the 
right). The reactor is initially empty; rates are scaled accord¬ 
ing to the volume-parameter Q = ~ 602 (see main 

text), which is the number of molecules at the fixed point, 
for both models and for both species. In the zero-dehciency 
case, stochastic dynamics only adds structure-less noise to the 
deterministic behavior. Instead, in the deficient case, while 
the deterministic system has damped oscillations towards the 
fixed point, oscillations are sustained in the corresponding 
stochastic dynamics, yielding a structural deviation between 
the two. 


to 


(Jt = 'E '^pi^t)Gp -\nxf (18) 

P 

As regards its stochastic counterpart, plugging the mass- 
action rates, Eq. 0, into Eq. ( p^ we obtain 

= J2(^p)tGp -J2^n\pt{X)X\]Lpt{X). (19) 

p X 

This is the first main result in our paper. Its most 
remarkable feature is that in the first term, related to 
the entropy flow to the environmenfP, only the “macro¬ 
scopic” average reaction velocity appears, and that “mi¬ 
croscopic” dependencies on X are within the second 
term, which is related to the system’s entropy change. At 
the trajectory level, this grants the validity of so-called 
Fluctuation Theorem^^, hence at is a proper notion of 
EPR. It is important, and a priori not obvious that the 
thermodynamic force Gp is the same in the stochastic 
and in the deterministic settings. 

Second, we define the correlation EPR as 5at := at —at 
and notice that, in the steady regime, it can be ex¬ 
pressed as a weighted difference between the average and 
the deterministic reaction velocity, as was anticipated in 
Eq. §■ Explicitly, we obtain a formula for the steady 
correlation EPR as a weighted sum of population mo¬ 


ments: 

5r^ oo 


[(^p) 


>(^oo) 


Go 


( 20 ) 


^ Gpkp ({X...{X-Up + l))oo - x'^£) .(21) 


The latter expression might pave the way for approxi¬ 
mate estimations of the correlation EPR based on Van 
Kampen’s system size expansion, moment-closure tech¬ 
niques or other diffusion approximations, provided due 
care is paid to the fact that such approximations often 
fail to reproduce the stochastic thermodynamics out of 
equilibriurrPH or even the distibution moment^^. 

Third, we evaluate the stochastic EPR when the sys¬ 
tem is in a product-form Poisson-like ensembl^with a 
generic time-dependent parameter yt^ 

Pois,,,(X) = —(22) 

with Zxo the normalization factor over Zxq • In this case 
it can be shown with few manipulations (see Appendix 
l^for a step-by-step derivation) that {vp )pois«, = Vp{yt), 
and consequently 


0-pois«, = 53 - Inyt • 53 ^pVpiyt)- (23) 

p p 


Notice that this expression coincides with the deter¬ 
ministic EPR at t ^ oc if the Chemical Master Equation 
admits a steady product-form Poissonian with parameter 
2/oo being a deterministic fixed point, and at all times if 
the system admits a product-form Poissonian with time- 
dependent parameter solving the deterministic rate equa¬ 
tions. 

Fourth, we investigate under which conditions such hy¬ 
pothesis are met. The ACK theorenP^ entails that, un¬ 
der our reversibility assumption, if the network has null 
deficiency, then the Chemical Master Equation admits a 
product-form Poissonian with parameter Xoo being the 
fixed point of the corresponding deterministic dynamics, 
which by Feinberg’s result^^ for ^ = 0 is unique and lo¬ 
cally stable. Hence the steady correlation EPR vanishes 
for zero-deficiency networks. For sake of reference we 
sketch a proof of the theorem in Appendix Further¬ 
more, it is known that in linear networks where no more 
than one molecule is consumed or produced at a time 
(i.e. Upp = 0,1), provided the system is prepared in 
a product-form Poissonian, it maintains such form at all 
times, with its parameter subjected to the correspond¬ 
ing rate equationPSl. Hence for linear CNs prepared in a 
product-form Poissonian ensemble, the correlation EPR 
vanishes at all times. These results thus generalize those 
by Mou et al.^^, who observed that the correlation EPR 
vanishes at steady states in linear networks. 


^ Notice that, because the range of summation is the lattice or- 
thant 2^Xo and not \X\ being the number of species, a 

“product-form Poisson-like” distribution is Poissonian in form 
but not in fact. 
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B. Numerical 

We will now illustrate the consequences of the ACK 
theorem and our findings with the aid of the above class 
of toy models. In fact we will further simplify the sce¬ 
nario by eliminating the waste W, which does not play 
any substantial kinetic role. Details on the simulation 
methods can be found in Appendix [C) 

Let be a scaling parameter regulating the system’s 
size and let x = N/Q be the concentration of N and 
y = EjQ that of E. A convenient choice of parameters 
is kp = ^ where K± are independent of the 

reaction, in their respective units (which depend on p). 
Then for given D all models turn out to have the same 
fixed point concentrations and steady EPR, making them 
easily comparable. Concentrations obey the system of 
rate equations 

x = K.- K_x - K+xy^ + K_y^^^ 

(24) 

y = n {K+xy^ - K_y^+^ + K_ - K^y^) . 

A fixed point is found at x^o = Voo = 1? for all values of 
m,n. Its stability depends on K-. The deter¬ 

ministic EPR at the fixed point is given by 

(Too = 3n{K+ - K_) In ^ (25) 

K — 

(notice that parameter Vi cancels within the logarithms, 
so that the EPR is extensive) and again it is independent 
of m, n. 

We will consider the cases n = 2, for values m = 
0,1, 2, 3, m = 0 being the zero-deficiency case, all others 
having 6=1. We take = 10, K- = 1 , which signi¬ 
fies that the system is very far from a detailed balanced 
thermodynamic equilibrium. We start from an empty 
reactor, xq = yo = 0. Eor these values the above fixed 
point is stable for all m < 4. Eor m = 0 the dynamics 
converges uniformly to the fixed point, as shown in the 
left-hand side of Eig.[2 A more interesting behavior ap¬ 
pears for higher m\ for m = 3 the deterministic system 
displays damped oscillations towards the fixed point (as 
shown by the innermost smoother lines in the left-hand 
side of Eig.[T). Indeed, for m = 4 the fixed point becomes 
unstable and the system displays steady oscillations. 

As regards the stochastic setting, so far our framework 
was that of ensemble thermodynamics, describing a large 
sample of processes at a given time. Erom now on we con¬ 
sider one given process in a large time. Indeed, Stochastic 
Thermodynamics has two complementary formulations: 
one along ensembles, and one along individual processed. 
The two frameworks are compatible, since the ergodic 
principle ensures that long-time averages almost surely 
(a.s.) equal ensemble averages at the steady state. In 
particular it can be proven that for the reaction velocity 

{vp)oo = lim a.s. (26) 

where is the number of times reaction p has been 

performed along the stochastic trajectory up to time t. 
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FIG. 2. The ACK theorem states that, if a CN has zero 
deficiency, then given an initial state (in our case, iV = E = 
0), the steady ensemble of the Chemical Master Equation has 
product form. Here we display color-plots of the histograms 
of the steady distribution of nutrient and energy molecules, 
for our toy models with m = 0,1, 2, 3, and rates scaled down 
by the volume-parameter D = 10~^^Na = 6.02, giving a low 
number of molecules at the steady state. Zebra lines (present, 
but not displayed for m > 0 for sake of better visualization) 
indicate that the stochastic dynamics preserves the parity of 
the energy molecules, which are produced in pairs. Owing to 
the outer smudge, the deficient models m = 1, 2, 3 have a non¬ 
product form distribution. The product-form distribution of 
the zero-deficiency case m = 0 is shown in more detail in 
Fig.ll 


Similarly, a histogram for the steady ensemble Poo (^: E) 
can be obtained by calculating the average time spent 
by the trajectory at state N, E. Let us then illustrate 
the ACK theorem. In Eig.[^ we provide color-plots for 
Poo{N,E). Eor m = 0, the color plot renders the distri¬ 
bution’s product-form. Zebra-lines are due to the fact 
that energy tokens are produced in pairs, hence starting 
from Xq = yo = 0 only even numbers of energy molecules 
can be populated. The same zebra-structure occurs for 
higher m > 0 , but for sake of better visualization we 
drew pixels twice the width, covering the whole area. 
The smudge in the color plots in Eig.j^for m > 0 reveals 
that the steady ensemble does not have product form. 
Instead, in the zero-deficiency case, Eig.[^ compares the 
histograms of the marginals for the energy and the nutri¬ 
ent, showing that they perfectly agree with the prediction 
from the product-form Poissonian. 

In Eig.j^we plot the average stochastic EPR as a func¬ 
tion of volume Cl. The perfect overlap between the de¬ 
terministic EPR (upper line) and the dots corresponding 
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Particle Numbers 


FIG. 3. The nonlinear CN 0 ^ N ^ 2E ^ 0 (corresponding 
to m = 0) has zero dehciency. Hence, by the ACK theorem its 
corresponding Chemical Master Equation affords a product- 
form steady ensemble, and the marginals for the number of 
nutrients and of energy molecules also have Poisson-like dis¬ 
tributions. We plot histograms for the populations of nutrient 
and energy molecules generated by stochastic simulations via 
Gillespie’s algorithm, with rates scaled by a volume parame¬ 
ter Q = 6.02, showing perfect agreement with the predictions 
of the ACK theorem. 


to the m = 0 case confirms our result that for deficiency- 
zero systems the correlation EPR vanishes. For m > 0 
this particular class of models has negative correlation 
EPR. The plots of the relative error in the inset show 
that the effect vanishes at large system sizes where fluc¬ 
tuations become negligible. 

Finally, another interesting aspect to inquire is the de¬ 
pendency of the correlation EPR on the affinity A = 
31ogiF+/iF_, which determines the distance from de¬ 
tailed balance, i.e. from thermodynamic equilibrium. In 
particular, we are interested in the so-called linear regime 
where the affinity is small and stationary currents are ap¬ 
proximately linear in the affinity. Then 

(5aoo = (^ - i)A^ (27) 

with the deterministic linear response eoeffieienti = (2/3. 
The inset in Fig.shows that in a model with nonvanish¬ 
ing deficiency, in the linear regime the correlation EPR, 
relative to the deterministic linear regime approximation, 
does not vanish in the limit A ^ 0, which implies that 
the stochastic linear response coefficient £ differs from the 
deterministic one. 

Our result proves that having (5 = 0 is a sufficient 
condition for a vanishing correlation EPR. A prelimi¬ 
nary question is then whether it is also necessary. The 
answer is trivially negative. In fact, if rates are such 
that detailed balance holds, then both the stochastic, 
the deterministic, and hence the correlation EPRs van¬ 
ish. More generally, for the ACK theorem to hold it 
is sufficient that the more general condition of eomplex 



EIG. 4. The main result of our paper is that dissipation 
(EPR) in stochasic chemical dynamics only coincides with 
the deterministic EPR when the GN has zero deficiency, and 
that already in simple systems intrinsic noise affects dissipa¬ 
tion. In the main frame we plot in log-log scale the stochastic 
EPR for our toy models, for all values m = 0,1, 2, 3, as a func¬ 
tion of the volume-parameter (2 that sets the average number 
of molecules present in the reactor at the steady ensemble. 
The up per straight line represents the deterministic value, 
Eq. (25). The dots on top of it are the values of the corre¬ 
sponding stochastic zero-deficiency system, m = 0. Models 
m > 1 with deficiency (5=1 have lower EPR than the deter¬ 
ministic model. An explanation for this is in Eig.[^ In the 
inset, we show that the relative error between stochastic and 
deterministic values decreases with volume. 


halanee holds: even if deficiency is greater than zero, 
rates can conjure in such a way that currents look “as 
if” the system had null deficiency. Furthermore, by the 
theory of Schnakenber^^ it can be shown that the cor¬ 
relation EPR can be decomposed in fundamental cycles 
(5cr(oc) = 5];^ [{Ja)oo - Ja{xoo)] with index (a span¬ 
ning a basis of the null space of the stoichiometric ma¬ 
trix, Aa a eyele affinity and a eyele eurrent. Cycle 
affinities are invariant under a wide range of transforma¬ 
tions of the rate constants which affect the cycle currents; 
hence even for non-complex balanced rates it might be 
feasible to tune the rates in such a way that several cycle 
contributions all cancel each other. 

The above argument rests on the fact that rate con¬ 
stants might be fine-tuned. The question becomes more 
interesting if properly reformulated. For systems with 
nonvanishing deficiency, complex-balanced rates are a set 
of measure zero in the space of possible rates. So, is 
the condition (5 = 0 necessary for a vanishing correlation 
EPR, for all possible values of rates? Very special sys¬ 
tems with nonvashing deficiency which still have Poisso- 
nian steady states have been foun(J^. An example is the 
chemical network X-j-Y 2X+Y, X 2X. In this case, 
the number of molecules of Y is constant and determines 
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FIG. 5. The affinity A = 31ogiF+/i^_ determines the dis¬ 
tance from thermodynamic equilibrium (detailed balance). In 
this figure we show the dependency of the deterministic and 
the stochastic EPRs with respect to the affinity, for m = 3 
and Q = 6.02, at fixed K- = 1 and variable iF+. The dashed 
curve is the linear-regime approximation of the deterministic 
EPR, where the current is approximately linear in the affin¬ 
ity and the EPR is approximated by a quadratic. Clearly the 
EPRs approach zero for vanishing affinity (no dissipation). 
The inset shows the error between stochastic and determin¬ 
istic EPR, relative to the linear approximation. The relative 
error increases with the affinity and, remarkably, it does not 
tend to vanish for A ^ 0. This implies that for nonvanishing 
deficiency, the Onsager coefficients of the deterministic and 
stochastic systems differ. 


the stoichiometric compatibility class where the dynam¬ 
ics is restricted. The deficiency is = 1, still the steady 
ensemble is a product-form Poissonian with parameter 
given by the solution of the deterministic equations of 
motion, and the correlation EPR can be easily shown 
to vanish. To take this class of cases into the descrip¬ 
tion, Cappelletti and Wiuf have introduced the concept 
of “stochastically complex-balanced” chemical reaction 
networks. The analysis of whether correlation EPR van¬ 
ishes for all values of the rates if and only if the network 
is stochastically complex-balanced goes beyond the scope 
of the present paper. 


IV. DISCUSSION AND CONCLUSIONS 

While it could have been expected that fluctuations 
would increase dissipation, our simple model displays the 
opposite behavior. This can be explained as follows. 
Notice that for m = 3 in Fig.[^ the stochastic dynam¬ 
ics has amplified oscillations, such as those characterized 
in Ref.l^, where a purely stochastic mechanism for bio¬ 
chemical oscillations was proposed. Such oscillations are 
forcedly stabilized in the deterministic setting. Hence 
the stochastic model is more flexible and capable of ex¬ 


ploring modes that the deterministic system abandons. 
Lower EPR then occurs when such modes are entropi- 
cally convenient. A way to characterize these modes is 
by a switching mechanism of chemical pathways. Fig.[^ 
details that in deficient networks, at low molecule num¬ 
bers certain reactions can be effectively shut off because 
of the temporary absence of a sufficient number of reac¬ 
tants. This phenomenon eventually reshapes the struc¬ 
ture of the irreversible closed reaction pathways that the 
system can locally perform. In our particular model, for 
low molecule numbers reaction 2 is inhibited, and the 
other two reactions alone do not contribute to dissipa¬ 
tion. Instead, in the CN with (5 = 0 the dissipative cycle 
can be performed at any particle number. 

The above example might then lead to hypothesize 
that the correlation EPR could be non-positive in gen¬ 
eral. This is not the case though. A counterexample can 
be found in the literature. The Schlogl model 0 ^ X, 
2X ^ 3X has deficiency (5 = 1, and its most important 
feature is that for certain critical values of the parameters 
it displays a bifurcation. Gaspard coimared stochastic 
and deterministic EPRs for this modeP^, and as can be 
observed from Fig. 2 in Ref.^^, close to the critical point 
the stochastic EPR is larger than the deterministic one, 
while in the bistable region it interpolates between the 
two possible values that the deterministic EPR takes at 
each of the two stable fixed points. 

Despite the fact that our toy model is oversimplified, 
the mechanisms we observed might carry out to more 
realistic networks. At the level of gene expression, it is 
known that intrinsic noise is a crucial factor in phenotypic 
variation within isogenic populations^. One step below, 
while in cells metabolites might be large in number, gene- 
expressed regulatory molecules might be very fevj^ allow¬ 
ing the switching mechanisms that we described above. 
In metabolism, the action of enzymes typically adds a 
level of complexity. In fact, most (if not all) of the reac¬ 
tions in biochemical CNs are not elementary, hence their 
connectivity and kinetic rules have to be determined a 
posteriori by advanced experimental methods (se^^ for 
a systematic review). Nevertheless, in our models the in¬ 
built deficient cycle could be seen as the core structure 
of any metabolic model. The network should be enriched 
by resolving individual metabolites within nutrients and 
waste, adding intermediate reactants such as cofactors 
and enzymes, resolving the environment and outer ther¬ 
modynamic cycles, separating time-scales and resorting 
to effective rate laws when applicable. As a proof of con¬ 
cept, all these operations will in general maintain the core 
cycle and hence the deficient character of the network, 
hence it can be argued that, because of its autocatalytic 
character, metabolism is deficient. 


^ In E. coli^ the lowest-concentration metabolite, nucleoside adeno¬ 
sine, is present in ~ 10^ copies, but over 80% of the variety 
proteins is much lower in copy numberJ^ 








To conclude, we emphasize that understanding ther¬ 
modynamic constraints on the regulation of met abolic 
networks is a crucial problem in CN reconstruct iorP^^. 
In this work we displayed a close connection between the 
topological notion of deficiency of a CN and nonequilib¬ 
rium thermodynamics, proving that at steady states only 
in zero-deficiency CNs the EPR evaluated by the mean- 
field deterministic theory coincides with that of the cor¬ 
responding stochastic model, accounting for stochastic 
variability in molecules’ number at low concentrations. 
For deficient CNs a nonvanishing correlation EPR quanti¬ 
fies the disagreement between deterministic and stochas¬ 
tic modeling, and at low molecule numbers this disagree¬ 
ment can be understood in terms of a switching mech¬ 
anisms of reaction pathways. A more detailed study of 
the conditions for positive vs. negative correlation EPR 
is demanded to future inquiry. Immediate perspectives 
also include the study of non-well-stirred mixtures, where 
react ion-diffusion processes allow for pattern formation, 
and of systems with separation of time scales and effec¬ 
tive enzymatic reactions. On the computational side, the 
more demanding stochastic techniques can be blended 
with deterministic algorithms to provide efficient tools 
for the systematic computation of the entropic balance 
of a CN, e.g. in software like COPASP^. More work 
has to be done to delineate future application of defi¬ 
ciency theory and stochastic thermodynamics to realistic 
metabolic networks. 
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Appendix A: Explicit derivation of Eq. ^ 23 ) 


From Eq. (14), plugging into the rates Eq. 0 and the 
Poisson-form distribution Eq. (22) we obtain 
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We now shift the summation over X to obtain 
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which is the desired result, Eq. (23). 
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FIG. 6. Chemical stochastic kinetics occurs on lattice or- 
thants, called stoichiometric compatibility classes (SCC). For 
our class of models, given an initial state, random jumps 
preserve the parity of the energy molecules (even or odd), 
hence there are two distinct SCCs. In the zero-deficiency 
case, m = 0, all of the drawn transitions are possible, and 
both SCCs can be obtained by repeatedly copy-pasting a mo¬ 
tif corresponding to the full CN, marked bold in the figure, 
through the whole lattice orthant. That is, locally each SCC 
looks like the full CN. Only cycling trajectories that carry 
a th ermod ynamic affinity contribute to the steady stochastic 
Epfpmii 

Hence, for m = 0, even for very low molecule num¬ 
bers it is always possible to perform the entropy-producing 
cycle. On the other hand, the structure of the SCCs for defi¬ 
cient networks is: for m = 1 dotted transitions type I are not 
feasible (since at least one energy token is needed to perform 
reaction p = -h2 and three energy tokens are needed to per¬ 
form p — —2), for m = 2 dotted transitions type I and II are 
switched off, and for m = 3 transitions type I, II and III are 
shut. Hence for low-enough molecule numbers, the stochas¬ 
tic trajectory explores a portion of the SCC where there is 
no possibility of producing entropy along an irreversible cycle 
(cycles consisting only of reactions p — ±1,±2 don’t dissi¬ 
pate). This explains the lower stochastic EPR observed for 
m = 1, 2, 3 in Fig.[^ 


Appendix B: Sketch of derivation of the deficiency-zero 
theorem 


One of the corollaries that incarnate the Anderson- 
Craciun-Kurtz theorenJ^ states that if a (weakly) re¬ 
versible CN has deficiency zero, then on each stoichiomet¬ 
ric compatibility classes the Chemical Master Equation 
admits a product-form Poisson-like steady distribution 
with parameter given by the unique fixed point of the 
corresponding rate equations. For sake of completeness, 
we provide the sketch of a derivation based on the graph- 
theoretical perspective that was briefly introduced in the 
main text. For another derivation based on quantum 
techniques, se^^. 
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Plugging the product-form Eq. (22) with parameter 
given by the fixed point cCoo into the generator Eq. (12), 
and using rates Eq. (11) one obtains 


LFois^^{X) = 
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(Bl) 


where we used Vp = i^-p — and antisymmetrized. We 
now observe that the sum over reaction vectors p > 0 can 
be commuted with a sum over complexes Y^, followed by 
a sum over all reactions p that have as a source com¬ 
plex. The latter information is stored into the incidence 
matrix d of the graph of complexes. Noticing that Up 
only depends on the complex of reactants ahead of p, we 
can write 


LPoiSa;^(X) = 


E 


Xp 


X- 




^ ^ ^^—piXoo) "^ApiXoo) • (B 2 ) 


After Eq. ( p!^ , the fixed point satisfies 

^ ^ ^p ['^-l-p(^oo) p(^oo)] — 0 


(B3) 


p>0 


which implies that u+p(cCoo) ~ is a right-null 

vector of the stoichiometric matrix. But if (5 = 0, then 
u+p(Xoo) — '^-p{xoo) is also a right-null vec tor of the in¬ 
cidence matrix (see last paragraph in Sec. IIB), hence 
Eq. (|B2|) vanishes. 


Appendix C: Materials and methods 


We employed the CN simulation software COPASP^to 
simulate the Chemical Master Equation via Gillespie’s al¬ 
gorithm, and the LSODA algorithm implemented in the 
scientific python stack (SciPy]P^ to solve deterministic 
rate equations. Histograms in Eig.j^and Eig.j^were sam¬ 
pled from stochastic trajectories for random-time change 
Markov jump processes spanning over 10^ s with a time 
resolution of 10“^ s, for a total of 10^ binned particle 
number pairs, while the stochastic time-courses in Eig|^ 
employ trajectories of 5 s with a resolution of 10“^ s. 
Each value for the average stochastic EPR in Eig.j^was 
calculated along single simulations of 10^ s. Notice that 
Gillespie’s algorithm keeps track of all reaction events, 
hence the final result for the stochastic average EPR is 
independent of time resolution. Eor the deterministic 
transients we used the same time-span and resolution as 
for the stochastic ones. The deterministic EPR was cal¬ 


culated via Eq. (25) and not from the simulation data. 


Thus it is only valid at the fixed point. 
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